// AWG_F
// Jared Bauters, June 2012
// University of California Santa Barbara

#include <mask.h>
#include <math.h>

// constants
#define PI 3.1415926535897932

// function
double free_propagation_region(double Ra,double alpha,double beta,double dR);
double awg_alcatel_const_wg_spacing_1_x_N_reverse(void);;

// Write AWG
int main(void) {
	
    setfont((char *)"../maskv6_00/caps.plf");
	msk_header("", "awg");
    
    awg_alcatel_const_wg_spacing_1_x_N_reverse();
    
	msk_trailer();
}

double free_propagation_region(double Ra,double alpha,double beta,double dR) {

    double aa,bb,cc,gamma,chi,ww;
    double R = Ra + dR;
    
    // Geometry
    gamma = PI - alpha - beta;
    aa = ((R)/2.0)*(sin(alpha)/sin(gamma));
    bb = ((R)/2.0)*(sin(beta)/sin(gamma));
    cc = sqrt( ((R)/2.0-aa)*((R)/2.0-aa) + ((R)-bb)*((R)-bb) - 2*((R)/2.0-aa)*((R)-bb)*cos(gamma) );
    chi = asin( (R-bb)/cc*sin(gamma) );
    ww = ((R)/2.0-aa);
    
    // Wedge of radius Ra and angle alpha
    push_cp();
    rotate(alpha);
    skip((R)/2.0);
    rotate(-PI/2);
    cw((R)/2.0,-2*alpha,(R));
    pop_cp();
    
    // Wedge of radius Ra/2 and angle beta
    push_cp();
    skip((R)/2.0);
    rotate(PI);
    rotate(beta);
    skip((R)/4.0);
    rotate(-PI/2);
    cw((R)/4.0,-2*beta,(R)/2.0);
    pop_cp();
    
    // Box 1
    push_cp();
    rotate(PI/2.0);
    pen_up();
    cw((R)/2.0,-beta,1e-6);
    pen_down();
    rotate(-(PI/2.0-chi));
    rect(cc,ww,1,3);
    pop_cp();
    
    // Box 2
    push_cp();
    rotate(-PI/2.0);
    pen_up();
    cw((R)/2.0,beta,1e-6);
    pen_down();
    rotate((PI/2.0-chi));
    rect(cc,ww,7,9);
    pop_cp();
    
    return RAD(chi);

}

double awg_alcatel_const_wg_spacing_1_x_N_reverse(void) {
    
    double x_axis = cp.a;
    
    // WG parameters
    double r_min = 1.5e-3;
    double wg_spacing = 15.0e-6;

    // AWG parameters
    int Nch = 8;
    int Na = 61;
    double Ra = 434.271e-6;
    double DL = 53.609e-6;
    double w = 2.8e-6;
    double wr = 5.0e-6;
    double wa = 5.0e-6;
    double dr = 1.8*wr;
    double da = wa + 0.4e-6;
    double dtheta = da/Ra;
    double io_l = 0.2e-3;
    double taper_l = 0.15e-3;
    double alpha = RAD(40.0);
    double beta = RAD(30.0);

    double theta_io = 1.64;
    double theta_min = theta_io-dtheta*(Na-1)/2;
    double r_0 = 1.815e-3;
    double s1_0 = 191.0e-6;
    double s2_0 = 0.0e-6;
    double P = 2*r_0*sin(theta_min/2);
    double awg_w = 2*( (Ra+s1_0)*cos(theta_min) + r_0*sin(theta_min) + s2_0 );

    // simplifying parameters
    double a,b,d,e,f,g,h,i;

    // arrayed waveguide iteration parameters
    double r_i,s1_i,s2_i,l_i,r_i_1,s1_i_1,s2_i_1,l_i_1,theta_i_1;
    
    layer(L_waveguides);
    
    pen_up();
    cw(r_min,-theta_io,w);
    taper(taper_l,w,wr);
    pen_down();
    
    // Receiver Waveguides
    push_cp();
    skip(Ra/2);
    rotate(PI);
    double dtheta_r = 2*dr/Ra;
    double theta_r_max = cp.a-x_axis - ((Nch/2.0-1.0))*dtheta_r - dtheta_r/2.0;
    double r_x_max = (Ra/2+taper_l)*cos(theta_r_max-PI)+r_min*sin(theta_r_max-PI);
    double r_p,r_x,r_aa,t1,t2;
    r_aa = cp.a;
    rotate((Nch/2.0)*dtheta_r + dtheta_r/2.0);
    for (int jj=0; jj < Nch; jj++) {
        rotate(-dtheta_r);
        
        r_x = r_x_max - (Ra/2+taper_l)*cos(cp.a-PI)-r_min*sin(cp.a+x_axis-PI);
        
        push_cp();
        skip(Ra/2);
        
        t1 = cp.a - r_aa;
        t2 = atan( sin(t1)/(cos(t1)+1) );
        rotate(t1-t2);
        
        push_cp();
        rotate(PI);
        sw(300e-9,wr);
        pop_cp();
        
        taper(taper_l,wr,w);
        fprintf(stderr,"x = %f um, y = %f um, a = %f rad \n",1e6*cp.x,1e6*cp.y,cp.a);
        //cw(r_min,(cp.a-PI),w);
        //sw(r_x,w);
        pop_cp();
    }
    pop_cp();
    
    // Input FPR
    free_propagation_region(Ra,alpha,beta,10e-9);
    
    // Arrayed Waveguides
    a = Ra;
    h = DL/2.0;
    i = wg_spacing;
    push_cp();
    rotate(dtheta*((Na-1)/2.0+1.0));
    for (int ii=0; ii < Na; ii++) {
        rotate(-dtheta);
        push_cp();
        skip(Ra);
        if (ii == 0) {
            b = cp.a-x_axis;
            r_i = r_0;
            s1_i = s1_0;
            s2_i = s2_0;
            l_i = r_min*b;
        }
        else {
            theta_i_1 = cp.a-x_axis-dtheta;
            b = cp.a-x_axis;
            d = s1_i_1;
            e = theta_i_1;
            f = r_i_1;
            g = s2_i_1;
            
            r_i = ( -1/(2-2*cos(b)-b*sin(b)) )*( (-1+cos(b))*(f+i-f*cos(e)-a*sin(b)+a*sin(e)+d*sin(e)) + sin(b)*(d+e*f+h+a*cos(b)-(a+d)*cos(e)-f*sin(e)) );
            s1_i = (a-d+b*f-e*f-h+b*i+(-a+d+e*f+h)*cos(b)-(a+d)*cos(b-e)+a*cos(e)+d*cos(e)-b*f*cos(e)-a*b*sin(b)-f*sin(b)-i*sin(b)+f*sin(b-e)+a*b*sin(e)+b*d*sin(e)+f*sin(e)) / (-2+2*cos(b)+b*sin(b));
            s2_i = d + f*e + g + h - s1_i - r_i*b;
            l_i = r_i*b + s1_i + s2_i;
        }
        //fprintf(stderr,"angle = %f\n r = %f mm, s1 = %f um, s2 = %f um, DL = %f um\n",b*180/PI,1e3*r_i,1e6*s1_i,1e6*s2_i,2e6*(l_i-l_i_1));
        
        taper(taper_l,wa,w);
        sw(s1_i-taper_l,w);
        cw(r_i,b,w);
        sw(2*s2_i,w);
        cw(r_i,b,w);
        sw(s1_i-taper_l,w);
        taper(taper_l,w,wa);
        r_i_1 = r_i;
        s1_i_1 = s1_i;
        s2_i_1 = s2_i;
        l_i_1 = l_i;
        pop_cp();
    }
    pop_cp();
    
    rotate(theta_io);
    skip(awg_w);
    rotate(-(PI-theta_io));
    
    // Output FPR
    free_propagation_region(Ra,RAD(40.0),RAD(30.0),10e-9);
    
    // Input Waveguide
    push_cp();
    sw(20e-9,wr); //overlap
    pop_cp();
    rotate(PI);
    taper(taper_l,wr,w);
    fprintf(stderr,"x = %f um, y = %f um, a = %f rad \n",1e6*cp.x,1e6*cp.y,cp.a);
    //cw(r_min,-theta_io,w);
    //sw(io_l,w);
    
    return awg_w;
}